// Copyright (c) 2009-2010, Jeremy Brewer
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.

package galaxie500.numeric;

import galaxie500.numeric.Function;

/**
 * Class implementing an adaptive integration algorithm
 * 
 * This class uses uses an adaptive integration strategy where the total
 * interval to evaluate is subdivided recursively. For each interval, it
 * evaluates 2 estimates of the integral over that interval, one using N points
 * and the other 2N points. If they differ by more than a threshold, the
 * interval is subdivided and the process repeated. Simpson's Rule is used to
 * evaluate the integral within each interval, and the overlapping function
 * calls between the N and 2N estimates are reused.
 * 
 * This class is based on the code and discussion from Algorithms in C++ by
 * Robert Sedgewick found in the integration chapter.
 * 
 * The integration routines take as input a class which implements the Function
 * interface. For simple functions, anonymous inner classes are very useful for
 * defining these classes.
 * 
 * Example Usage:
 * 
 * <pre>
 * Function gaussian = new Function() {
 *   &#064;Override
 *   public double valueAt(double x) {
 *     return (1.0 / Math.sqrt(2.0 * Math.PI)) * Math.exp(-0.5 * x * x);
 *   }
 * };
 * Integrator integrator = new Integrator();
 * double result = integrator.integrate(gaussian, -1.0e6, 1.0e6);
 * </pre>
 * 
 * @author Jeremy Brewer
 */
public class Integrator {

  /** Error tolerance (absolute) used as convergence criterion. */
  private double maxError = 1.0e-6;

  /** The number of intervals (N) used each level of recursive integration. */
  private int numIntervals = 5;

  /** The number of recursive splits accumulated during integration. */
  private int numSplits = 0;

  /** The maximum allowed number of splits allowed. */
  private int maxNumSplits = 1500;

  /**
   * Temporary variable to hold the integral estimate with numIntervals points
   */
  private double integral = 0.0;

  /**
   * Temporary variable to hold the integral estimate with 2 * numIntervals
   * points
   */
  private double integral2 = 0.0;

  public Integrator() {
  }

  // Accessors.

  public double getMaxError() {
    return maxError;
  }

  /**
   * Controls tolerance for convergence. Experimentation seems to indicate that
   * a simple absolute error check works best (why I'm not sure).
   * 
   * @param maxError Maximum allowed tolerance (absolute)
   */
  public void setMaxError(double maxError) {
    assert maxError > 0.0 : "Invalid value for maxError: " + maxError;
    this.maxError = maxError;
  }

  public int getNumIntervals() {
    return numIntervals;
  }

  /**
   * At each level, we evaluate the integral using numIntervals and 2 *
   * numIntervals to decide when to recurse further. For Simpson's rule, there
   * are 4 * numIntervals + 1 function calls per level. Thus, setting this
   * determines the minimum number of functions calls to evaluate an integral.
   * Sometimes, for a very wide interval, you might not sample enough points to
   * get a good answer (say for a very narrowly peaked function). In those
   * cases, changing this value to something higher can help, though in general
   * you don't need to alter the default value.
   */
  public void setNumIntervals(int numIntervals) {
    assert numIntervals > 1 : "Invalid numIntervals: " + numIntervals;
    this.numIntervals = numIntervals;
  }

  public int getNumSplits() {
    return numSplits;
  }

  /**
   * Sets the maximum number of recursive splits allowed in the integration.
   * 
   * @param maxNumSplits The maximum number of recursive splits allowed
   */
  public void setMaxNumSplits(int maxNumSplits) {
    assert maxNumSplits > 0 : "Invalid value of maxNumSplits: " +
                              maxNumSplits;
    this.maxNumSplits = maxNumSplits;
  }

  public int getMaxNumSplits() {
    return maxNumSplits;
  }

  // Metadata methods.

  /**
   * Returns the number of function calls per interval.
   * 
   * @return Number of function calls per recursive interval
   */
  public int callsPerLevel() {
    return 4 * numIntervals + 1;
  }

  /**
   * This method can be invoked after integrating to determine how many function
   * calls were required.
   * 
   * @return The total number of function calls used in an integration.
   */
  public int numFunctionCalls() {
    if (numSplits == -1) return 0;
    return 2 * numSplits * callsPerLevel() + callsPerLevel();
  }

  /**
   * Adaptive integration routine. Evaluates the integral of function f from a
   * to b using an adaptive algorithm which places more points where f varies
   * the fastest for increased efficiency.
   * 
   * @param f The function to integrate
   * @param a The starting point of integration
   * @param b The ending point of integration
   * @return The value of the integral of f from a to b
   */
  public double integrate(Function f, double a, double b) {
    numSplits = -1; // Special case for equal endpoints.
    if (a == b) return 0.0;
    numSplits = 0;
    return integrateImpl(f, a, b);
  }

  /**
   * Implementation of the adaptive integration routine which assumes a != b.
   * 
   * @param f The function to integrate
   * @param a The starting point of integration
   * @param b The ending point of integration
   * @return The value of the integral of f from a to b
   */
  private double integrateImpl(Function f, double a, double b) {
    integral = 0.0;
    integral2 = 0.0;
    simpsonIntegrate(f, a, b, numIntervals);

    // This is a simple absolute error check. I experimented with using
    // relative floating point comparisons as well as scaling the error by the
    // current interval size, but this approach seems to work best (better error
    // in fewer function calls). We should probably do more thorough testing to
    // determine the optimum splitting strategy here.
    if (Math.abs(integral - integral2) > maxError) {
      double mid = 0.5 * (a + b);
      numSplits++;
      if (numSplits > maxNumSplits) {
        throw new RuntimeException("More than " + maxNumSplits + " splits");
      }
      return integrateImpl(f, a, mid) + integrateImpl(f, mid, b);
    }

    return integral2;
  }

  /**
   * Internal integration routine. Integrates using Simpson's rule and stores
   * the estimate for the integral of f from a to b for both numIntervals and 2
   * * numIntervals in private variables integral and integral2 respectively.
   * The function calls between both estimates are shared.
   * 
   * @param f The function to integrate
   * @param a The starting point of integration
   * @param b The ending point of integration
   * @param numIntervals The number of intervals to use for integration
   */
  private void simpsonIntegrate(Function f, double a, double b,
                                int numIntervals) {
    // Width of each subinterval.
    double dx = (b - a) / (double) (numIntervals);
    double dxHalf = 0.5 * dx;
    double dxQuart = 0.25 * dx;
    double dxThreeQuart = 0.75 * dx;

    // Contribution from the first interval.
    double sum = 0.0;
    double sumHalf = f.valueAt(a + dxHalf);
    double sumQuart = f.valueAt(a + dxQuart) + f.valueAt(a + dxThreeQuart);

    // Contribution from the inner intervals.
    for (int i = 1; i < numIntervals; i++) {
      double x = a + i * dx;
      double xHalf = x + dxHalf;
      double xQuart = x + dxQuart;
      double xThreeQuart = x + dxThreeQuart;
      sum += f.valueAt(x);
      sumHalf += f.valueAt(xHalf);
      sumQuart += f.valueAt(xQuart) + f.valueAt(xThreeQuart);
    }

    // Contribution from the endpoints.
    double fa = f.valueAt(a);
    double fb = f.valueAt(b);
    integral = (dx / 6.0) * (fa + fb + 2.0 * sum + 4.0 * sumHalf);
    integral2 =
        (dx / 12.0) * (fa + fb + 2.0 * (sum + sumHalf) + 4.0 * sumQuart);
  }
}
